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Introduction 


The  dynamic  compaction  of  loose  dry  granular  material  is  fundamentally  a  multi-scale  problem. 
Many  engineering  problems  involving  rapid  1  oading  conditions  exist  at  the  bulk  scale,  and  exam  pies 
include  planetary  impact  and  crater  formation,  tectonic  plate  movement,  ballistic  im  pact  and 
penetration  events.  On  a  bulk  scale,  which  includes  thousands  of  grains  of  sand,  the  material  behaves 
somewhat  like  a  homogeneous  hydrodynamic  continuum.  However  as  the  characteristic  dimension  of 
the  analysis  is  reduced,  the  do  minant  phenomenological  behavior  is  altered.  At  smaller  scales,  on  the 
order  to  tens  of  grains  of  sand,  the  material  form  s  heterogeneous  structure,  such  as  stress  bridgin  g, 
which  serves  to  distribute  loads.  The  result  is  a  co  mplicated  stress  field  where  grains  within  a  bridge 
network  carry  nearly  all  of  the  load  and  grains  outside  of  a  bridge  network  carry  nearly  no  load.  At  yet 
smaller  scale,  on  the  order  of  the  grain  contact  surface,  the  material  contact  tensor  determines  the 
transmittal  of  stress  from  one  grain  to  the  next.  At  even  smaller  scale  (those  less  than  the  grain  itself) 
the  non-isotropic  material  rheolog  y  behavior,  such  as  grain  y  ield  and  defo  rmation,  plasticity  and 
twinning,  establish  the  dominant  phenomenological  behavior  characteristics.  Phenomenology  at  all  of 
these  scales  feeds  back  to  the  bulk  scale  and  ultimately  dictates  the  response  of  the  overall  bulk  system. 
However  it  remains  unclear  as  to  whic  h  mechanisms  dominate  and  which,  if  any,  might  therefore  be 
neglected.  The  si  mulations  presented  here  resolve  t  he  dynamic  compaction  behavior  from  the  bulk 
scale  (although  not  relative  to  most  engineering  systems)  to  the  grain,  but  not  sub-grain.  At  the  scale  of 
the  grain  and  smaller  the  material  is  treated  as  a  horn  ogenous  continuum  solid.  Thus  t  he  term  meso 
scale,  which  means  middle  scale,  is  used  to  describe  the  simulations  presented  here. 

In  addition  to  dynamics  at  multiple  scales,  the  variability  associated  with  loading  rate  can  affect  the 
dynamic  response  of  materials.  Quartz  has  famously  demonstrated  a  rich  variety  of  behaviors  as  a  result 
of  shock  loading,  including  phase  transitions,  mat  erial  property  variability  as  a  function  of  driving 
pressure  and  sam  pie  thickness,  as  well  as  stress  r  elaxation  effects  [i-iv].  The  strain  rates  under 
investigation  here(<103  s'1)  are  modest  compared  to  shock  loading  ( 1 0  5~106  s'1)  but  clearly  in  the 
dynamic  regime.  These  loading  rates  are  more  characteristic  of  conditions  found  further  way  from  an 
impact  event,  where  the  bulk  material  does  not  necessarily  experience  uniform  loading  in  excess  of  the 
Hugoniot  elastic  limit  (HEL).  An  exam  pie  of  such  a  loading  condition  is  at  ransmitted  stress  wave 
emanating  from  a  land  mine  explosion.  The  analysis  is  complicated  by  the  heterogeneous  nature  of  the 
material  and  therefore  the  loading.  Although  the  bulk  stress  state  of  the  sand  does  not  exceed  the  HEL, 
this  does  not  imply  that  i  ndividual  grains  won’t  exceed  the  HE  L.  The  response  of  heterogeneous 
systems  when  subjected  to  these  in  termediate  strain  rates  r  emains  of  fundam  ental  importance  in 
completing  our  understanding  of  the  dynamic  behavior  of  sand. 

Thus,  our  aim  is  to  explore  the  dy  namic  compaction  of  a  heterogeneous  sy  stem  using  a  mesoscale 
numeric  approach  at  low  to  m  oderate  strain  ra  tes.  Ultim  ately  it  is  hoped  that  these  mesoscale 
simulations  can  be  used  to  augmen  t  experimentation  by  exploring  the  state  space  of  a  heter  ogeneous 
system  were  experimental  measurements  are  difficult  or  i  mpossible,  and  serve  as  a  platform  by  which 
continuum  models  can  be  developed  or  refined. 
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Background 


Due  to  the  computational  expense  associated  with  resolving  a  mesoscale  simulation,  the  usefulness 
as  applied  to  large  sc  ale  engineering  systems  remains  limited.  Thus  the  need  for  continu  um  modeling 
persists.  The  variation  in  materi  al  properties  associated  with  heterogeneous  systems,  the  inclusion  of 
void  space  and  void  collapse,  coupled  with  the  lack  of  adequate  experimental  techniques  to  probe  t  he 
heterogeneous  state,  make  continuum  modeling  challenging.  Usually,  granular  materials  are  modeled 
as  a  continuum  (i.e.  the  grains  and  porosity  are  treated  as  a  n  ew  single  material)  and  a  re  assigned 
average  properties.  In  so  doing  the  heterogeneous  nature  of  the  material  and  the  grain  interactions  are 
lost.  In  these  continuum  approaches  additional  constitutive  relations  such  as  the  sphere  collapse  models 
[v],  snowplow  [vi],  P-a  [mi],  P-A  [niu]  or  27-alpha  [ix]  models  are  used  to  govern  the  removal  of 
porosity. 

Developments  in  computer  architecture  and  soph  istication  in  par  allelization  have  allowed  large- 
scale  simulations  to  explore  multi  scale  pheno  mena.  A  host  o  f  researchers,  utilizing  a  variety  of 
numeric  formulations  such  as  Eulerian  h  ydrocodes,  finite  element  or  discrete  element  methods,  have 
successfully  used  mesoscale  simulations  to  m  odel  many  types  of  high  rat  e  dynamic  phenomena 
involving  heterogeneous  materials,  including  the  shock  compaction  of  elemental  metals  and  alloys,  the 
behavior  of  energetic  materials,  ceramic  alloys  and  earth  materials  such  as  sand  [x,xi,xii].  These  studies 
have  demonstrated  that  mesoscale  simulations  can  be  used  to  s  uccessfully  describe  the  high  rate 
compaction  dynamics  without  the  use  of  additional  constitutive  relations,  thus  corroborating  predictions 
from  various  analy  tic  porous  collapse  models  and  experimental  observations  [  xiii,xiv,xv].  In  most  of 
these  studies  the  focus  was  to  num  erically  predict  the  bulk  material  behavior  at  high  strain  rates,  i.e. 
Hugoniot  behavior.  It  was  shown  that  si  mple  formulations  and  simple  descriptions  on  a  grain  level  can 
combine  to  resolve  complicated  dynamics  observed  on  the  bulk  scale  such  as  the  formation  of  dynamic 
force  chains,  void  collapse,  turbulent  like  behavior  and  hot  spot  formation.  This  study  differs  from 
these  in  that  we  aim  to  explore  the  dynamic  compaction  at  much  lower  strain  rates,  near  103  s'1. 


Experimental  Data 


The  experimental  data,  which  is  compared  to  the 
simulations  presented  he  re,  was  obtained  from  a 
split-Hopkinson  pressure  bar  [xvi-xix].  The  material 
used  in  these  experimental  investigations  was  silica 
based,  kiln  dried  fine  grai  n  sand.  The  particle  size 
distribution  was  between  200-450  pm.  The  dynamic 
compressive  response  of  the  sand  was  investigated 
at  various  moisture  contents  rangin  g  from  3%  to 
20%  by  weight  with  all  specimens  having  a  dr  y 
density  of  1.50  g/cm3.  The  specimens  were  confined 
using  a  hard  ened  4340  s  teel  tube  with  an  outer 
diameter  of  25.4  mm,  inner  diameter  of  1  9.1  mm 
and  length  o  f  50.8  mm .  The  steel  tube  is  used  to 
achieve  high  confinement  pressures  and  to  replicate 
a  uni-axial  st  rain  condition.  The  incident  bar  was 


Strain 

Figure  1:  Uni-axial  stress-strain  behavior  of 
sand,  as  well  as  the  incident  bar  velocity, 
obtained  from  a  split  Hopkinson  bar  at  a  strain 
rate  of  606  s-1  [16-18]. 
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well  controlled  in  order  to  produce  a  near  linear  strain  (constant  strain-rate  of  approximately  600  s’1) 
over  the  test  interval  of  0.2  milliseconds.  A  representative  example  of  the  resulting  axial  stress  data  i  s 
presented  in  Figure  1  as  a  function  of  strain.  The  reduction  in  stress  at  a  strai  n  of  approximately  0.2 
represents  significant  unloading  of  the  sample.  The  strain  profile  was  differentiated  in  order  to  produce 
an  incident  bar  velocity  profile,  also  illustrated  in  Fig.  1.  This  velocity  profile  was  used  as  a  boundary 
condition  for  the  simulations  presented  below. 


Numeric  Simulations 

The  simulations  presented  here  were  performed  using  the  CTFl  hydrocode  [xx].  The  sand  was 
modeled  either  as  a  collection  of  quartz  spheres  in  a  three-dimensional  rectilinear  dom  ain  for  t  he 
mesoscale  simulations  or  as  a  single  representative  material  in  a  one-dimensional  domain  for  the 
continuum  simulations.  This  was  done  so  that  the  performance  of  each  type  of  simulation  could  be 
assessed  and  compared.  For  both  types  of  simulations,  the  motion  of  a  rigid  driver  plate  re  presenting 
the  incident  bar  was  prescribed  on  the  left  side  of  the  domain  and  a  symmetry  boundary  condition  was 
imposed  on  the  right  side  of  the  domain  to  represent  the  longitudinal  center  of  the  sam  pie.  Figure  1 
presents  the  velocity  versus  strain  history  that  was  used  to  prescribe  the  motion  of  the  incident  bar  for 
these  simulations.  For  the  three-dimensional  simulations,  periodic  boundary  conditions  were  i  mposed 
in  the  lateral  direction.  Greater  detail  regarding  each  type  of  simulation  is  presented  below. 

One-Dimensional  Continuum  Calculations 

Of  fundamental  importance  when  conducting  continuum  simulations  of  heterogeneous  materials  is 
selecting  the  right  constitutive  relations;  in  this  wor  k  several  were  investigated,  however  onl  y  one  is 
presented  here:  a  tabular  EOS  (sesame  table)  with  a  P-a  compaction  model  [vii,xxi,xxii].  This  equation 
of  state  had  been  previously  developed  for  simulations  of  buried  land  mine  explosions.  The  amplitude 
and  duration  of  trans  mitted  stress  waves  far  fro  m  land  mine  explosions  is  si  milar  to  t  he  loading 
conditions  applied  by  the  split  Flopkinson  bar.  A  su  ite  of  sesame  tables  for  vary  ing  soil  moisture 
content  has  been  included  in  the  stan  dard  release  of  CTFl  since  version  8.0.  For  the  simulations 
presented  here  the  eos7860  tabular  sesame  was  used,  which  was  formulated  for  4%  moisture  and  bulk 
Si02  sand,  where  the  grain  densit  y  is  2.65  g/cc  and  the  bul  k  density  is  1  .56  g/cc.  The  underlying 
reference  curve  for  the  construction  of  these  tables  was  the  shock  Flugoniot  for  quartz  measured  by 
Wackerle  [i].  Kerley  suggested  that  the  pore  compaction  pressure  be  approximate  with: 

P  =  (500  MPa)e  192w,  (1) 

where  w  is  the  moisture  content  of  the  sand  [xxii  ].  Thus  assuming  4%  m  oisture,  we  obtain  a  pore 
compaction  pressure  of  232  MPa,  however  it  was  found  that  a  value  of  1 ,440  MPa  best  fit  the 
Flopkinson  bar  data.  The  bulk  strength  of  the  sand  was  modeled  with  a  pressure  dependent  geological 
yield  surface  (GEO),  previously  implemented  in  CTFl  [xxiii].  The  formulation  for  this  yield  surface  is 
presented  in  Equation  2  and  is  illustrated  in  Fig.  2.  Experimental  data  from  buried  land  mine  explosions 
suggest  the  yield-pressure  correlation,  dY/dP,  is  approxim  ately  2  and  would  remain  positive  through 
compaction.  The  yield  offset  strength,  Y,  is  approximately  80  MPa,  and  the  bulk,  zero  pressure,  yield 
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strength,  Y0,  should  be  between  100-300  kPa,  [  xxi,xxii,xxiv],  The  remaining  material  properties  , 
obtained  from  a  vari  ety  of  sources  [xxv-xxvii],  are  listed  in  Table  1.  Note  that  variations  in  the  bulk 
yield  strength  had  little  effect  on  the  axial  stress  obtained  given  the  loads  far  exceed  these  strengths. 


Pressure,  P  (MPa) 

Figure  2:  Pressure  dependent  geological  yield 
surface 
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Table  1:  Baseline  constitutive  relation  constants  for  the  continuum  simulations 


Tabular  Equation  of  State 

Parameters'1 

Value 

Pressure  Dependent  Geological  Yield 
Parameters'1 

Value 

Bulk  Density,  p00  [g/cc] 

1.56 

Max.  Strength,  Y  [MPa] 

300 

Grain  Density,  p  [g/cc ] 

2.56 

Yield  Surface  Slope,  dY/dP 

2 

Pore  compaction  pressure,  Ps  [MPa] 

232 

Yield  strength  at  zero  pressure,  Y0  [MPa] 

0.3 

Zero  stress  shock  speed,  c0  [km/s] 

5.99 

Poisson  ratio,  v 

0.32 

Hugoniot  slope,  s 

2.345 

Tensile  strength,  cts  [GPa] 

0.0001 

Kerley  [xxi,xxii] 


Three-Dimensional  Mesoscale  Calculations 

The  computational  domain  and  m  aterial  properties  utilized  for  the  m  esoscale  simulations  are 
presented  in  Fig.  3  and  T  able  2  respectively.  Unlike  the  one-dimensional  co  ntinuum  simulations,  a 
mesoscale  approach  resolves  the  grain  interaction  directly,  thus  there  is  no  need  to  approximate  the  bulk 
material  properties  or  include  a  pore  compaction  model. 

A  Mie-Griineisen  equation  of  state  (EOS)  [xxviii]  with  a  perfectly  elastic -plastic  strength  m  odel 
was  used  to  m  odel  each  quartz  sand  grain.  The  EO  S  parameters  are  based  on  cry  stalline  quartz,  i.e. 
silicon  dioxide  (Si02);  the  main  constituent  of  Ottawa  sand  is  9  9.8%  pure  quartz  [  xxix].  Quartz  has 
famously  demonstrated  a  rich  variety  of  behaviors  as  a  result  of  shock  loading,  includ  ing  phase 
transitions,  material  property  variability  as  a  function  of  driving  pressure  and  sample  thickness,  as  well 
as  stress  relaxation  effects.  Given  the  nature  of  the  Eulerian  calculations,  i.e.  the  inability  to  specify 
crystallographic  orientation  of  indi  vidual  sand  grains  as  well  as  the  polycrystalline  nature  of  grains, 
different  bulk  material  properties  are  assigned  to  eac  h  sand  grain.  For  example,  Table  2  list  the  shock 
Hugoniot  properties  of  quartz  as  a  function  of  t  he  orientation  of  the  crystalline  structure  [ii].  However, 
instead  of  assigning  averaged  properties  to  all  of  the  grains,  we  explore  the  effect  of  assigning  the  x-cut 
shock  Hugoniot  properties  to  half  of  th  e  grains  in  the  simulation  and  z-cut  properties  to  the  other  half. 
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The  effect  of  vary  ing  the  material  properties  is  not  new  to  this  work;  a  more  detailed  exploration  of 
theses  effects  can  be  found  in  Branson,  Wells  and  Strack  [  xxx]  and  Borg  [xxxi],  The  specific  heat  of 
quartz  is  reported  between  0.74  J/g/K  and  0.85  J/g/K  [xxxii,xxxiii]. 

For  these  simulations,  the  grains  behave  elastically  until  the  specified  dynamic  strength  is  achieved, 
at  which  point  the  material  plastically  flows.  Estimates  of  the  Hugoniot  elastic  limit  (HEL)  for  the  x- 
and  y-cut  axes  are  between  5.5  GPa  to  8.5  GPa  and  10  to  15  GPa  for  the  z-cut  axis,  however  hydrostatic 
data  indicate  the  lower  limit  is  closer  to  2.3  GPa  [i,iii].  Simulations  were  conducted  in  which  the  HEL 
for  all  the  grains  was  varied  from  2  to  12.4  GPa.  The  measured  Poisson’s  ratio  for  qu  artz  varies 
between  0.11-0.15,  however  it  was  set  to  0.15  for  these  simulations  [i]. 

The  measured  dynamic  tensile  strength,  i.e.  spall  strength,  of  quartz  exhibits  variations  with  respect 
to  loading  conditions.  Kanel  et.  al.  report  that  the  spall  strength  for  x-cut  quartz  reaches  4  GPa  for 
loadings  below  the  HEL  and  falls  to  zero  for  1  oading  conditions  near  the  HEL  [xxxiv] .  The  d  ynamic 
tensile  strength  is  extremely  high  as  compared  to  the  static  fracture  of  44  MPa  [xxxv].  The  loss  in  spall 
strength  is  thought  to  result  from  accumulated  damage,  i.e.  cracking  of  brittle  single  cry  stals  under 
compression,  ultimately  resulting  in  a  total  loss  of  te  nsile  strength  near  the  HEL.  Much  like  the  shock 
Hugoniot  parameters,  and  in  order  to  understand  the  effect  of  strength  on  the  bulk  compaction,  a  range 
of  spall  strengths  were  assigned  to  the  grains  within  a  given  simulation.  Simulations  were  carried  out 
where  the  tensile  strength  was  varied  from  pre-damaged  or  static  fracture  strengths,  near  zero,  to  4  GPa. 


Figure  3:  Three  dimensional  view  of 
mesoscale  geometry.  Variations  in  gray  scale 
are  for  visualization  purposes  only 


Table  2:  Baseline  material  and  constitutive  constants  for 
the  mesoscale  simulations 


Parameter 

Quartz 

Density,  p  [ g/cm 

2.65 

Zero  stress  shock  speed,  C0  [km/s\ 

x-cut 

5.610 

z-cut 

6.329 

Hugoniot  slope,  s 

x-cut 

1.07 

z-cut 

1.56 

Griineisen  coefficient,  T =V(dP/dE)v 

0.9 

Specific  heat,  Cv  [J/(g-K)\ 

0.85 

Bulk  dynamic  yield  strength,  Y  [GPa\ 

x-cut  (low,  average,  high) 

4.1,  5.8,  7.0 

z-cut  (low,  average,  high) 

8.2,  10.3,  12.4 

Poisson’s  ratio,  v 

0.15 

Fracture  strength,  <ys  [GPa 

0.044-  15  GPa 

Given  experimental  observations  that  friction,  be  it  grain-on-grain  or  between  the  sand  and  test 
fixture,  has  a  significant  effect  on  the  resulting  stress-strain  behavior,  we  included  some  form  of  friction 
in  the  m  esoscale  simulation.  Given  the  Eulerian  nature  of  CTH,  im  posing  a  co  mplete  contact  stress 
tensor  on  the  surface  of  each  grain  would  be  difficult  at  best.  H  owever,  CTH  does  allow  t  he  user  to 
specify  how  the  strength  of  a  mixed  material  cells  is  calculated.  When  setting  up  the  com  putational 
domains,  care  was  taken  to  assign  adjacent  grains  a  different  material  number,  thus  enabling  some 
control  over  the  grain  contact  strength  in  mixed  (i.e.  touching)  cells.  We  explored  the  effects  of  mixed 
cells  being  assigned  either  zero  or  the  volum  e  fraction  averaged  bulk  yield  strength:  sliding  or  stiction 
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respectively.  Sliding  m  eans  there  is  no  strength  between  materials  in  m  ixed  cells,  where  as  stiction 
means  that  the  initial  grain  network  is  welded.  In  either  case,  once  the  stress  in  a  given  cell  exceeds  the 
assigned  strength,  the  material  flows  at  the  specified  yield  stress. 

In  previous  high  strain  rate  studies,  10  6  s’1,  conducted  by  the  authors  it  has  been  found  t  hat  mesh 
convergence  is  achieved  for  10  computational  cells  per  grain  per  spatial  direction.  In  this  st  udy,  where 
the  strain  rates  are  much  lower,  the  mesh  was  such  that  20  cells  per  grain  were  needed.  Given  these 
requirements  and  the  three-dimensional  nature  of  the  domain,  the  simulations  required  approximately 
100  hours  utilizing  64  pr  ocessors,  whereas  the  one -dimensional  simulations  required  15  minutes 
utilizing  a  single  processor.  One  can  see  the  value  of  further  developing  continuum  level  simulations. 


Results 

Mesoscale  Simulation  Results 

In  order  to  e  xplore  the  mesoscale  results,  the  longitudinal  stress,  crxx,  was  averaged  in  the  lateral 
direction  for  a  given  longitudinal  position.  The  results  are  presented  in  Fig.  4  at  two  different  snapshots 
in  time,  which  translates  to  19  and  31%  strain,  for  both  the  stiction  and  sliding  grain  conditions.  The 
simulations  indicate  that  sliding  results  in  substantia  lly  softer  bulk  stress  as  compared  to  stiction.  A 
particular  feature  of  these  relatively  slow  strain  rate  events  is  that  the  stress  state  is  fairly  constant  ahead 
of  the  incident  bar  (driver  plate).  The  transmittal  of  information  occurs  at  a  rate,  on  the  order  of  the 
sound  speed  of  quartz  (6,  000  m/s)  or  possibly  as  slow  as  the  sound  speed  o  f  bulk  sand  ( —100  m/s), 
either  of  which  is  much  faster  than  that  of  the  dr  iver  plate  (~1  m/s).  In  the  context  of  a  hy  drocode 
written  specifically  to  resolve  wave  phenom  ena  this  translates  into  hundreds  of  stress  r  everberations 
between  the  sample  thickness.  Thus  as  the  driver  plate  co  mpacts  the  s  and  at  relatively  constant 
velocity,  a  time  varying  (strain  varying)  stress-strain  state  is  achieved  within  the  sand. 

These  results  provide  a  justification  to  average  the  entire  stress  state  of  the  sand  in  order  to  obtain  a 
bulk  axial  stress  state  for  a  given  displacement,  i.e.  strain.  The  results,  presented  in  Fig.  5,  indicate  that 
the  sand  with  sliding  grain  contact  requires  very  little  applied  stress  to  strain  to  nearly  25%  as  compared 
to  the  grains  with  stiction.  The  simulations  indi  cate  that  the  bulk  com  paction  snow  plows  in  its 
response,  i.e.  significant  pore  removal  for  little  to  no  applied  stress,  before  the  bulk  response  stiffens.  It 
is  interesting  to  note  that  the  maximum  pack  density  of  a  monodisperse  collection  of  spheres  is  n/  VTs  or 
74%.  Thus  grain  packing  without  friction  slides  t  ogether  rather  unimpeded  until  it  is  mechanically 
locks  at  near  25%  strain,  at  which  point  the  response  stiffens.  This  behavior  is  illustrated  in  Fig.  5. 

There  are  two  modes  in  which  grain  material  can  be  rearrange:  1)  centroid  motion,  i.e.  the  grains 
moving  without  deforming  and  2)  grain  deformation,  i.e.  a  grain  being  deformed  into  a  g  ap.  For  both 
the  sliding  and  stiction  simulations  m  aterial  is  rearranged  through  a  combination  processes  1)  and  2). 
The  stiction  simulations  appear  to  have  less  c  entroid  motion  and  more  grain  deformation  as  compared 
to  the  sliding  simulation.  The  sliding  simulations  do  not  achieve  a  face  centered  cubic  arrangement,  i.e. 
maximum  pack  density.  These  assessments  are  purely  qualitative;  it  is  difficult  to  quantify  how  much 
of  the  material  rearrangement  is  due  to  deformation  and  how  much  is  due  to  centroid  rearrangement. 


6 


Figure  4:  Average  longitudinal  stress  for  grains  with 
either  stiction  or  sliding  at  two  different  strains 


Figure  5  Stress-strain  behavior  for  various  grain 
and  strength  configurations 


When  compared  to  experimental  stre  ss-strain  data,  both  the  stiction  and  slidi  ng  grain  conditions 
under  predict  the  experimental  behavior.  Thus  the  average  dynamic  yield  strength  was  increased  by  an 
order  of  magnitude,  from  44  MPa  to  440  MPa,  to  ass  ess  the  effects  on  the  bulk  response;  the  results  are 
presented  in  Fig.  6  along  with  experimental  results  for  comparison.  The  results  indicate  an  increase  in 
the  dynamic  yield  strength,  over  that  of  the  st  atic  yield  strength,  is  necessary  to  repr  oduce  the 
experimental  results.  Since  the  sliding  grain  co  nditions  indicated  near  sn  owplow  like  behavior, 
whereas  the  data  did  not,  sliding  simulations  with  higher  yield  strengths  are  not  presented. 


Strain 

Figure  6:  Bulk  stress-strain  for  the  mesoscale 
simulations  as  a  function  of  grain  boundary  and 
dynamic  yield  strength 


Strain 

Figure  7:  Stress  tensor  components  for  the  mesoscale 
simulations  with  a  dynamic  yield  strength  of  440  MPa 
and  stiction  grain 


In  order  to  gain  a  better  appreciat  ion  for  th  e  complete  stress  stat  e  achieved  within  the 
heterogeneous  system  both  the  averaged  lateral  str  ess  and  the  averaged  absolute  value  of  the  shear 
stress  were  computed  from  the  mesoscale  simulations  and  the  results  are  presented  in  Fig.  7.  Since  the 
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average  shear  stress  is  zero;  the  absolute  value  of  the  shear  stress  was  computed  to  gain  an  appreciation 
for  amplitude  of  the  relative  shear  compared  to  the  longitudinal  and  lateral  stress. 


Lateral  Stress,  Ow  (MPa)  Shear  Stress,  CTX v  (MPa) 


(a)  (b)  (c) 

Figure  8:  Stress  distribution  of  averaged  longitudinal,  lateral  and  shear  stress  at  0.05,  0.1  and  0.15  strain 


In  order  to  better  chara  cterize  the  state  of  stress  within  the  sand  mixture,  stress  distributions  were 
obtained  from  the  mesoscale  data.  Figure  7  presents  the  distribution  of  the  various  tensor  components 
of  stress  for  all  the  com  putational  cells  within  the  s  imulation,  at  three  different  strains.  The  resulting 
distributions  are  non-normal  (i.e.  non-Gaussian).  There  are  a  large  number  of  computational  cells  that 
experience  no  applied  stress;  thes  e  are  represented  as  a  large  frequency  of  occurrence  at  zero.  As  a 
result,  the  most  probable  occurrence  of  stress,  i.e.  the  average  str  ess,  does  not  appear  centered  on  the 
distribution.  For  example,  at  0.15  strain  Fi  g.  7  indicates  the  average  longitudinal  stress,  axx,  is 
approximately  160  MPa,  which  appear  s  lower  than  the  peak  occurrence  of  stress  presented  in  Fig.  8a 
for  a  strain  of  0.15.  In  add  ition  there  are  relatively  high  stress  hot  spots  experienced  within  the  sand 
which  are  not  completely  characterized  by  the  average,  i.e.  some  computational  cells  are  experiencing 
stress  nearly  an  order  of  magnitude  higher  than  th  e  average.  Insights  such  as  these,  which  are  obtained 
from  mesoscale  simulations,  could  provide  the  way  forward  in  our  understanding  of  rapidly  compacted 
heterogeneous  systems,  including  both  theoretical  development  and  suggested  experimental  validation. 

One-dimensional  Simulation  Results 


Figure  9  presents  the  stress  strain  results  from  the  one-dimensional  simulations  described  above. 
Like  the  three-dimensional  simulations,  the  one-dimensional  simulations  with  the  original  param  eters 
under  predict  the  dynamic  response  of  the  sand.  By  increasing  the  stress  at  which  all  of  the  pores  ar  e 
removed  (i.e.  the  stress  at  which  the  P-a  curve  intersects  the  quartz  Hugoniot)  from  440  MPa  to  1440 
MPa,  the  sim  ulations  more  accurately  follow  the  experimentally  measured  material  response.  This 
modification  of  the  baseline  one-dimensional  material  parameters  is  not  in  complete  agreement  with  the 
three-dimensional  mesoscale  simulations.  In  the  mesoscale  simulations  the  co  mplete  removal  of 
porosity  corresponds  to  t  he  intersection  of  the  bulk  dynamic  response  and  the  underl  ying  material 
Hugoniot,  which  is  near  the  dynam  ic  strength  assigned  to  the  underlying  material,  i.e.  440  MPa.  Thus 
the  one-dimensional  simulations  require  additional  considerations. 
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Figure  9:  Bulk  stress-strain  curve  for  the  continuum 
simulations  yield  strength  with  varying  pore 
compaction  strength 


Figure  10:  Stress  tensor  components  for  the 
continuum  simulations  with  a  pore  compaction 
strength  of  1500  MPa 


Discussion 

One  unanticipated  aspect  of  this  work  was  th  e  higher  resol  ution  required  to  co  nverge  the 
simulation,  which  is  nearly  twice  that  of  previous  higher  strain  rate  mesoscale  simulations  [xii],  In 
general  it  is  accepted  that  resolving  each  grain  with  10  computational  cells  per  grain,  per  axial  direction 
is  sufficient.  However  t  hese  lower  s  train  rate  si  mulations  demonstrated  significant  variations  in 
response  up  to  20  cells  per  grain  per  axial  direction.  This  could  be  a  result  of  the  high  num  her  of 
iterations  required  to  run  the  simulation  out  to  nearly  6  ms  and  the  advection  problems  as  sociated  with 
the  large  number  of  iteration  performed  in  Eulerian  simulations. 

For  the  simulations  presented  here,  a  distribution  of  yield  and  fracture  strength  was  applied  to  the 
grain  network  in  order  to  acce  ss  the  e  ffect  on  the  bulk  response.  When  appl  ying  a  distribution  of 
strength  the  key  consideration  is  the  average  value.  For  the  loading  confi  gurations  investigated  here, 
and  the  distribution  of  yield  strengths  from,  2  to  12.5  GPa,  it  was  found  that  if  the  average  did  not  vary 
then  the  bulk  response  did  not  var  y.  This  is  not  the  case  for  higher  strain  rat  e  investigations,  10 6  s"1 
[xxxi].  In  the  higher  strain  rate  studies  it  was  f  ound  that  an  un  der lying  skeletal  network  of  higher 
strength  grains  can  form  a  percolated  stress  bridge  which  can  stiffen  the  bulk  response,  this  is  especially 
true  for  stress  levels  near  full  compaction.  For  the  lower  strain  rate  studies  conducted  here,  no  such 
phenomena  occurred.  We  are  not  sure  if  this  is  due  to  the  allowable  time  for  rearrangement,  the  low 
stress  levels  induced  or  the  relatively  narrow  range  of  yield  strength  prescribed.  Simulations  presented 
here  were  sensitive  to  the  variations  in  fracture  strength,  however.  Varying  the  distribution  of  fracture 
strength  could  alter  the  bulk  response  of  the  grain  network,  even  though  the  average  value  remained  the 
same.  This  i  s  especially  true  if  one  ap  plies  a  dramatic  distribution  of  fracture  strengths,  i.e.  half  th  e 
grains  having  a  fracture  strength  of  zero  and  half  the  grains  having  a  fracture  strength  nea  r  the  HEL. 
These  results  require  further  investigation. 

The  data  presented  in  Fig.  7  is  re -plotted  in  Fig.  11a  long  with  two-dimensional  mesoscale  results 
from  high  strain  rate  sim  ulations  and  a  ccompanying  experiments,  shown  as  open  circles  a  nd  dashed 
lines  respectively  [x,xxxvi].  It  i  s  interesting  to  note  that  at  higher  strain  rates,  1 06  s"1,  the  presence  of 
stiction  over  predicts  the  stiffness  of  the  response  as  co  mpared  to  experim  ental  data  [  x,xxxvi] .  It  is 
important  to  keep  in  mind  that  the  high  strain  rate  experimental  data  was  obtained  from  a  one¬ 
dimensional  flyer  plate  experiment,  whereas  the  data  presented  above  was  obtained  from  a  Hopkinson 
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bar.  The  mesoscale  simulations  are  nearly  identical  except  they  are  at  higher  incident  bar  velocities  and 
two-dimensional. 


Figure  11:  A  comparison  of  high,  106  s-1,  and  low,  103  s-1,  strain  rate  data  and  simulation.  Low  strain 
rate  data  from  Fig.  7.  High  strain  rate  data,  open  circles,  and  accompanying  mesoscale  simulations,  dashed 

lines,  with  stiction  or  sliding  [10,36] 


As  Fig.  11  indicates,  the  sliding  boundary  condition  better  represents  the  ex  perimentally  determined 
dynamic  response  of  the  materi  al  at  h  igher  strain  rate.  Both  high  and  low  strain  rate  simulations 
suggest  that  grain  contact  interactions  can  have  a  significant  effect  on  the  bulk  response  of  the 
heterogeneous  material  but  that  the  relevant  contrib  ution  of  friction  varies  as  the  strain  rate  varies,  i.e. 
friction  is  important  at  low  strain  rates,  with  dimensioning  importance  at  higher  strain  rates.  Exploring 
these  effects  will  continue  for  future  studies. 


Conclusions 

Mesoscale  simulations,  such  as  those  presented  here,  provide  a  rich  backdrop  for  the  exploration  of 
rapidly  compacted  heterogeneous  systems.  Mesoscale  simulations,  which  are  a  fairly  new  idea  within 
the  computational  community,  have  yet  to  realize  their  full  potential  as  a  research  tool  in  conjunction 
with  experimental  techniques,  theoretical  development  and/or  continuum  simulations.  The  simulations 
presented  here  indicate  that  low  strain  rate  Eulerian  mesoscale  simulations  require  fairly  high  resolution 
to  accurately  predict  the  averaged  longitudinal  stress  of  a  dynamically  compacted  sand  sam  pie.  In 
addition,  the  simulations  required  that  t  he  dynamic  strength  be  an  order  of  m  agnitude  higher  than  the 
static  yield  strength.  The  simulations  further  indicate  that  the  presence  of  stiction  (friction)  is  necessary 
to  accurately  predict  the  stress-strain  behavior  as  compared  to  experimental  data.  From  the  simulations 
a  complete  description  of  the  averaged  stress  state,  the  longitudinal,  lateral  and  shear  stress,  can  be 
obtained.  In  addition,  the  simulations  also  indicate  that  this  stress  state  is  non-normal  distributed  within 
the  sand  sample  and  hot  spots  of  stress  can  exists  which  are  nearly  ten  times  that  of  the  av  erage  axial 
stress  measured  by  the  experiment. 
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